//
// Torino2009.cxx
//
// Created by TB on 11/23/11.
// Copyright 2011 BNL. All rights reserved.
//

#include "tppmc/sivers/Torino2009.h"

#include <cstdlib> // For abs(int)
#include <sstream>
#include <vector>

#include <boost/foreach.hpp>

#include "tppmc/Enumerations.h" // For LHAPDF::Flavour
#include "tppmc/Grv98Pdf.h"
#include "tppmc/Particle.h"
#include "tppmc/PartonContent.h"

namespace tppmc {

   namespace sivers {

      Torino2009::Torino2009(bool bound)
      : mUseBound(bound)
      , mUnpolarised(Grv98Pdf::LO) {
         mParameters.SetReadFile("tables/sivers_parameters.dat");
         mParameters.ReadParameters();
         mParameters.UseMainSet();
      }

      Torino2009::~Torino2009() {
      }

      // total(x, kperp) = unpol(x, kperp) + 0.5 * DeltaF(x, kperp) * spinfactor
      // |spinfactor| < 1
      void Torino2009::Evaluate(const pythpart& parton,
                                double scale,
                                PartonContent& partcontent) const {
         // Light quark flavours are supported:
         int flavours[] = {
            LHAPDF::UP, LHAPDF::DOWN, LHAPDF::STRANGE,
            LHAPDF::UBAR, LHAPDF::DBAR, LHAPDF::SBAR
         };
         // The Sivers extraction is parameterised in terms of the GRV98
         // unpolarised PDF set.
         // Evaluate the unpolarised distribution.
         PartonContent unpolarised(0.);
         mUnpolarised.Evaluate(parton, scale, unpolarised);
         
         partcontent = PartonContent(0.); // Initialise to zeros
         
         BOOST_FOREACH(int i, flavours) {
            partcontent[i] = unpolarised[i];
            // If using the positivity bound, we just return the unpolarised,
            // otherwise multiply by the fit parameters.
            if(not mUseBound) {
               // The parameterisation for DeltaF, which is added to the unpolarised
               // term, is 2 * N(x) * h(kperp) * unpolarised(x, kperp)
               // N(x) * h(kperp) is returned by fit->EvaluateSivers().
               // Note that the 2 above cancels out the 0.5 before the DeltaF term.
               // Don't actually have unpolarised(x, kperp) just unpolarised(x)
               const FitValues* fit = mParameters.parameters.Values(i);
               if(not fit) {
                  continue;
               } // if
                 // Fit values multiply the unpolarised distribution (already set)
               partcontent[i] *= fit->EvaluateSivers(parton);
            } // if
         } // BOOST_FOREACH
           // The d quark sivers shoudld be negative.
           // If using the parameterisation it is already negative from.
         if(mUseBound) {
            partcontent[LHAPDF::DOWN] = -partcontent[LHAPDF::DOWN];
         } // if
      }

      bool Torino2009::Supports(const pythpart& parton) const {
         return ::abs(parton.Pdg() < 4) or 21 == parton.Pdg();
      }

      std::string Torino2009::Name() const {
         std::string name("Anselmino et al. 2009 Sivers extraction");
         if(mUseBound) {
            name.append(" at positivity bound");
         } // if
         return name;
      }
   } // namespace sivers
} // namespace tppmc
